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Abstract 

We present an efficient quantum algorithm for beyond-Born-Oppenheimer molecular en¬ 
ergy computations. Our approach combines the quantum full configuration interaction method 
with the nuclear orbital plus molecular orbital (NOMO) method. We give the details of the 
algorithm and demonstrate its performance by classical simulations. Two isotopomers of the 
hydrogen molecule (H2, HT) were chosen as representative examples and calculations of the 
lowest rotationless vibrational transition energies were simulated. 
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1 Introduction 


Exact computations and simulations of quantum systems on a classical computer are computa¬ 
tionally hard. This stems from the fact that the dimensionality of the Hilbert space needed for 
the description of a studied quantum system scales exponentially with its size. One of the con¬ 
sequences is e.g. the prohibitive exponential scaling of the full configuration interaction (FCI) 
method. Quantum computers^, on the other hand, offer an exponential speed-up for this task 
as was first noticed by Feynman and Manin- ®. The underlying idea, which employes mapping of 
the Hilbert space of a studied quantum system onto the Hilbert space of a register of quantum bits 
(qubits), both of them being exponentially, large, can in fact be adopted also in quantum chemistry. 

The past few years have witnessed a remarkable interest in the application of quantum com¬ 
puting for solving of different problems in quantum chemistry. Among others, quantum algorithms 
for non-relativistic®® as well as relativistic — molecular FCI energy calculations, quantum chem¬ 
ical dynamic^®, or calculations of molecular properties 17 were developed. For a complete list of 
relevant papers, we refer the reader to recent reviews 18 20 . Efficient quantum chemical simulations 
are indeed believed to belong to the first practical applications of quantum computers. This is also 
supported by recent proof-of-principle few-qubit experiments 21 26 . Several improvements reducing 
the resource requirements of fault-tolerant implementation and thus paving the way for practical 
simulations were presented in 27 . 

In this paper, motivated by the fact that non-Born-Oppenheimer (non-BOA) effects play an 
essential role in wide range fields (e.g. proton tunnelling in DNA damage), we generalise the 
applicability of the quantum FCI algorithm 1113 for beyond-BOA computations. We achieve this 
by combining qFCI with the NOMO method 28 '" 34 . We should however note that our attempt is 
not the first one dealing with beyond-BOA computations on a quantum computer. In®, it was 
shown that simulating all electron-nuclear and inter-electronic interactions (and thus going beyond 
BOA) is somewhat surprisingly faster and more efficient than BOA for systems with more than 
four atoms. Nevertheless, the aforementioned approach is based on the first quantized formulation, 
thus completely different from ours. 

The structure of the paper is following. In Section [2] we shortly review the basic concepts 
of the phase estimation-based quantum FCI algorithm 11 ®^, in Section [ 3 ] we do the same for the 
NOMO method, and in Section [4] we present details of our quantum algorithm for beyond-BOA 
computations, which is a combination of both approaches. The performance of the proposed scheme 
is presented in Section [5] by classical simulations of H 2 and HT energy computations. 

2 Quantum FCI algorithm 

An efficient quantum FCI (qFCI) algorithm for calculations of nonrelativistic molecular energies 
employing the phase estimation algorithm (PEA) of Abrams and Lloyd' was proposed in the pio¬ 
neering work by Aspuru-Guzik et al . 11 . It was later simplified by replacing of PEA with its iterative 
version, iterative phase estimation algorithm (IPEA)®®. (I)PEA is a quantum algorithm for ob¬ 
taining the eigenvalue of a unitary operator U, based on a given initial guess of the corresponding 
eigenvector. Since a unitary U can be written as U = e iH , with H Hermitian, the (I)PEA can be 
viewed as a quantum substitute of the classical diagonalization. 

Suppose that \u) is an eigenvector of U and that it holds 

U\u)=e 2 ^\u), ^ 6(0,1), (1) 

where cf is the phase which is estimated by the algorithm. In case of the original PEA, a quantum 
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register is divided into two parts. The first one is the read-out part composed of m qubits on which 
the binary representation of the estimate of phase cp is eventually measured. After the application of 
Hadamard gates in the read-out part of the quantum register followed by application of a sequence 
of controlled-H 2 operations (fcs are nonnegative integers from 1 to m), the register is transformed 
into 

i«s> = 4 = m 

V 3=0 V j=0 

The next part of the algorithm is the inverse quantum Fourier transform (QFT) 1 performed on the 
read-out part of the register. The whole register is transformed into \2 m (p)\u) and the phase can 
be extracted from its first part. 

Iterative version, IPEA adopts the ideas of measurement-based quantum computing® and 
reduces the computational resources by using only single read-out qubit. If (p is expressed in the 
binary form: cp = 0 .</>i </>2 ... (p m , (pi = {0,1}, one bit of (p is measured on the read-out qubit at each 
iteration step. The algorithm is iterated backwards from the least significant bits of cp to the most 
significant ones. The fc-th iteration is shown in Figure [lj 



<Pk 


Figure 1: The k -th iteration of the iterative phase estimation algorithm (IPEA). H denotes 
Hadamard gate (7r/2 rotation) and the feedback angle Uk depends on the previously measured bits 
(see Eq. [Ij). 


The equivalent of QFT is a single qubit ^-rotation R z , whose angle ujk depends on the results of 
the previously measured bits 


Rz {f^k ) — 

1 0 \ 

^ 0 e 2niu}k ) 

(3) 

Wfc = 

m-k -\-1 , 

tyk+i-l 

Z-, 2 * ’ 

i=2 

(4) 


followed by a Hadamard gate. 

Depending on how the second part of the quantum register is treated in between individual 
iterations, we distinguish two versions of IPEA 11 . In case of version A, it is maintained during 
all iterations (initialised only once). The biggest advantage of this approach is that one always 
ends up with one of the eigenstates of U. On the other hand, the quantum coherence is needed for 
the whole algorithm which makes this version more difficult for an experimental realization. IPEA 
version B is on contrary characterised by reinitialization of the second part of the register at every 
iteration step. Therefore, the quantum coherence is needed only within each iteration separately. 

The (I)PEA algorithm can be exploited for ab initio quantum chemical calculations, if we take 
U in the form 7 11 


U = e irA , 


(5) 
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where H is the molecular electronic Hamiltonian (up to now, only Born-Oppenheimer Hamiltonians 
have been considered) and r is a suitable parameter which assures </> being in the interval (0,1). 
The electronic Hamiltonian can be expressed in the second quantized form a^® 


H = 


Zz 

pq 


L pq Uj p'' 


2 ^ 

pqrs 


( 6 ) 


where h pq and V pqrs are one and two-electron integrals in the molecular spin orbital basis and a\ 
and at are fermionic creation and annihilation operators. Since these operators in general do not 
commute, exponential of the Hamiltonian cannot be written as a product of exponentials of indi¬ 
vidual hx , but a numerical approximation must be used®. The first-order Trotter approximation 37 
has the form 


e irH = e iTj2 L x=1 hx = ^ JJ /hxr/N^ + Q^/ N y (7) 

X=1 

When representing the quantum chemical wave function on a quantum register, the simplest 
approach (but the least economical one in terms of number of qubits) is so called direct mapping. 
It directly assigns individual spin orbitals (or in relativistic case Kramers pair bispinors) to qubits, 
since they can be either occupied or unoccupied (occupation number basis), corresponding to |1) 
or |0) states. Jordan-Wigner transformation (JWT)®® is then used to express fermionic operators 
in terms of Pauli a matrices. JWT has the form 



where a± = 1/2 (a x ± ia y ) and the superscript denotes the qubit on which the matrix operates. 
Alternatively, the Bravyi-Kitaev transformation 39 40 which balances locality of occupation and 
parity information and reduces the simulation cost from 0(n) to O(logn) for one fermionic operator 
when compared to JWT may be used. We would like to note that more compact mappings e.g. 
from subspace of fixed-electron-number wave functions or spin-adapted wave functions can also be 
used efficiently in connection with quantum sparse simulation algorithms*- 41 ! 

Regarding the overall scaling of the qFCI algorithm, Wecker et a/! 42 found that the computa¬ 
tional time for bounded error scales with the number of spin orbitals N as 0(N 9 ) on average and 
as 0(N 11 ) at worst. Using the Bravyi-Kitaev transformation 39 40 instead of JWT would decrease 
the scaling to 0(N 8 log N) or O(lV 10 log./V) respectively. Poulin et a/J 43 used testing set of real 
molecules and observed even more feasible scaling of the Trotter-Suzuki time step leading to the 
overall scaling 0(IV 4 ~ 5 ' 5 ). The computational cost bounds can be further improved considering 
local molecular basis sets®! 

The qFCI algorithm 1113 requires an initial guess of the exact eigenstate, whose quality in¬ 
fluences the success probability of measuring the desired energy. This can be either a classical 
approximation [e.g. complete active space (CAS) based wave function 1 ®^, an exact state pre¬ 
pared by the adiabatic state preparation method*- 11 * 45 ! or by the algorithmic cooling method 4 -^, 
or also a unitary coupled cluster approximation optimised by the recently presented combined 
classical-quantum variational approach® 47 [ 

In order to increase the overall success probability of (I)PEA, the whole algorithm is repeated 
and the correct phase cj) decided from the majority voting. In case of version B, individual iterations 
are independently repeated and correct <pk values decided from the majority voting. 
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3 NOMO method 


The nuclear orbital plus molecular orbital (NOMO) theory^ 28 34 (other authors denote similar the¬ 
ories with different acronyms, e.g. ENMO 48 ^ or NEO^ is an extension of MO theory to the 
non-BOA problem, which employs the idea of a nuclear orbital (NO), as a one-particle orbital of 
a nucleus. Since electrons and nuclei are treated on equal footing in NOMO framework, it goes 
beyond the Born-Oppenheimer and adiabatic approximation^®. 

In the NOMO method, Gaussian-type functions are adopted for electronic as well as nuclear 
basis functions, which leads to difficulties in gauging the total-energy accuracy because of the poor 
description of translational and rotational motions. For this reason, translation- and rotation-free 
(TRF) approach has been developed by eliminating the contribution of translation and rotation 
from the total Hamiltonian 31 34 . The TRF Hamiltonian has the form 


#TRF = H - T T - Tr 

= ^TRF + ^TRF + ^TRF + ^TRF + ^TRF • (9) 

The total translational Hamiltonian Tt is simply subtracted, but since molecular rotations and 
vibrations are coupled, the contribution of rotation cannot be entirely eliminated. The only viable 
way how to subtract rotation is by Taylor expansion of the rotational Hamiltonian Tr with respect 
to A(as defined in®). In our numerical study (see Section [ 5 ]), this was done just to zeroth order. 

The Hamiltonian Htrf contains one-particle (nucleus: T FRF , electron: T FRF ) and two-particle 
(nucleus-nucleus: R frf , nucleus-electron: RF , electron-electron: R frf ) terms and has the same 
second-quantized form as in Eq. [6j If, for simplicity, we consider now only one kind of nucleus and 
use big subscripts {P, Q, . ..} for NOs and small {p, q, . ..} for MOs, we can write 


Htrf = ^2 h e p e q ala q + ^ h n p Q a) p aQ 

pq PQ 

+ 2 ^ r pqrs^pO’q^sa r + - ^2 VpQRsaV a Qasa R 

pqrs PQRS 

+ V pQrsaW Q a s a r . ( 10 ) 

pQrS 

In principle, the NOMO/FCI theory for a complete configuration space is an exact theory. 
In practice however, due to the exponential scaling of classical FCI, some approximation has to 
be adopted. Different kinds of NOMO post-Hartree-Fock methods analogous to those from the 
conventional Born-Oppenheimer electronic structure theory has been developed, to name a few e.g. 
NOMO/MP2 30 ®, NOMO/Ci 5 ®, or NOMO/CC® 1 . One of the drawbacks of the NOMO theory 
is a slow convergence of the n-e correlation effect with respect to CI/CC expansion 52 54 . As we 
are dealing with the NOMO/FCI theory here, which is, as already mentioned, an exact theory, we 
avoid this problem. For more details about the NOMO methodology and the discussion of its pros 
and cons, we refer the reader to the original literature 31 52 . 


4 Quantum NOMO/FCI algorithm 


In this section we elaborate how the qFCI algorithm is employed for the NOMO Hamiltonian (10). 
Despite restricting ourselves to distinguishable nuclei in the proof-of-principle numerical simula- 
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tions, the presentation in this section is intentionally kept as general as possible, i.e. considering 
several types of fermionic and/or bosonic nuclei. 

The NOMO Hamiltonian (10) can be cast to the general form Q, if the indices p,q,r,s are 
allowed to run over nuclei spin orbitals as well, i.e. h pq and V pqrs stands for one and two-particle 
integrals over molecular and/or nuclear spin orbitals. For this we define the following ordering of 
the electronic and nuclear spin orbitals. Let us consider a general system consisting of K kinds of 
different nuclei, where we have Nq molecular spin orbitals for no electrons, N\ nuclear spin orbitals 
for ni nuclei of the 1st kind, etc., up to Nk nuclear spin orbitals for nx nuclei of the A-th kind. 
Say that A/ erm kinds correspond to fermionic nuclei and the rest to bosonic ones. 

Therefore p , q, r , s E {1, 2,..., Nt}, where Nt is the total number of spin orbitals 


I< 

N T = J2^k. ( 11 ) 

k =0 

The operator a\, in (JH]) denotes the creation operator for p- th spin orbital [creating either electron 
(for p < No) or nucleus (for p > No)] and a q is the annihilation operator. 

The k- th kind nuclear spin orbital indices belong to the set Sy. 


Sk — { 4 ) Ik + 1 ; • ■ • 5 Ik + Nk — 1 }, 


where 4 is the k -th nuclei starting index 


( 12 ) 


fc-i 


4 = ^ Ni + 1 . 


(13) 


4 = 0 


All creation and annihilation operators commute with both creation and annihilation operators 
for different particles 


= = [“p> “<?] = °> ( 14 ) 

if p E Sk, q E Si and k / l. 

If particles of k-th kind are fermionic, the usual anticommutation relation holds 


— {dp,dq} — 0, {at, a,} — 5 pq , 
if p, q E Sk and fermions, 

while bosonic nuclei must obey the commutation relations 


(15) 


[aj,, at] — [a p , a q 


= 0, 
if p, q E S k 


[d.p, a q 


= 


•pqi 

and bosons. 


(16) 


Distinguishable identical nuclei can be considered as separate classes of different nuclei, each class 


consisting of a single nucleus, so that the relations (14) apply. 


The mapping between (anti)symmetrized products of spin orbitals of our general system and 
states of a quantum register can be constructed in the following way. Let us order the nuclei classes 
so that fermionic nuclei (k < Kf eTm ) precedes the bosonic ones (k > Kf eim ) and discuss the mapping 
for different particle types separately. 
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4.0.1 Fermions 

For fermions of &:-th kind (p E Sk), the standard Jordan-Wigner mapping ([8]) can be used, which 
in our index convention can be expressed as 



i<p 


(17) 


(18) 


When individual exp (ihxT/N) (see Eq. [7| terms are implemented, this approach leads to the 
computational cost O(N^). Alternatively, the more economic Bravyi-Kitaev transformation®®] 
[G(N% log Nk)] can be employed. 


4.0.2 Bosons 


In contrast to fermions, more than one qubit is needed to store an occupation number of a bosonic 
spin orbital. When considering bosonic particles of k- th kind, the occupation number f(p) of the 
spin orbital p can acquire values from 0 to n*,. The so called direct boson mapping® uses rq, + 1 
qubits to store this occupation number in the following way 


' /( p )- 1 


•n-k 


I fip)) = 


|0>i 


| 0 >, 


(19) 


® I !>/(*) ® 

\ t=o / Vi=/(p)+l / 

For clarity, we take into account only the part of a quantum register that corresponds to p-th spin 
orbital (in bosonic case one creation or annihilation operator does not act on qubits corresponding 
to different spin orbitals). 

From (16) follows the action of creation and annihilation operators 


«pl/(p)> = Vf(p) + + 1 ), 

(20) 

a P \f(p)) = Vf(p)\f(p) ~ !)> 

(21) 

a p 0) = 0. 

(22) 


As the maximum occupation /(p) is the number of particles n^, one more condition has to be 
introduced 


a{K> = 0. (23) 

In analogy to Eqs. and ([l8|) , Somma et al. 5 proposed the direct boson mapping of the form 


nk- 1 


^ Vj + 1 (o+®<^- + 1 J. 

(24) 

j= o 


n k -l 


^ V / 7 + 1 ( (T -®°'+ +i J- 

(25) 

3=0 
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Following their complexity analysis, we can express the overall computational cost of a single Trotter 
step of the general system as 


O 


E N 2 k N?\og(N k Ni) 


+ 


. 0<k<l<K te , 


+o 


+o 


E 


n k n l^k^l ] + 


, Kferm < Ck<l<K 


E log(^V fc )iV/A ^ fc 2 


(26) 


- 0<fc<A' ferm 

A'ferm <1<K 


where individual terms correspond to contributions from fermion-fermion, boson-boson and fermion- 
boson interactions, respectively, and N k < Ni for each 7Ff erm +1 < k < l is supposed for simplicity. 


The Bravyi-Kitaev transformation is considered for fermions in (26), as well as in the following 


formulae (30) and (32). 


In contrast to the direct boson mapping, we propose the compact boson mapping which uses 
only [dog 2 (nfc + 1)] qubits to represent the binary expansion of f(p), as shown below (the least 
significant bit is written leftmost) 


o' 

II 


|0)|0)|0).. 

-|0), 


1 f(p) = 1 ) 

i— y 

|1)|0)|0).. 

-|0), 


I/(P) = 2> 

i y 

|0)|1)|0).. 

-|0), 


s 

II 

^co 


|1)|1>|0).. 

,. 0),etc. 

(27) 


The representation of creation and annihilation operators for p -th spin orbital is defined implic¬ 
itly by relations (20)-(23) and we can express them semi-formally by the formulae 


n k ~ 1 

E Vi + i 1 f(p)=j +1) (f(p)=j\, 

7—0 

(28) 

n/g — 1 

E Vi + 1 1 f(p)=j) (f{p)=j + 1|- 

(29) 


3=0 


Rather than using formulae (28) and (29) explicitly, circuit representation of appropriate com¬ 
binations of the bosonic creation and annihilation operators occuring in (10) were investigated. In 
the Appendix we show that this approach leads to the overall computational cost of a single Trotter 
step 













where 


0[ J2 N%N?log(N k N l )\ + 

V 0<k<l<I< ielnl ) 

S y' j Ng (n k ,ni)N k N 2 j+ 

erm <k<l<K / 


+o( ]T Ag(l ; ni) \og{N k )NiNl j, 


■ CKAKifferm 

K{erm<l<K 




( 30 ) 


1 


1 


N g (n, m) = y + -xy{x + y)+ 

+|®V - ]-xy A + -*-y 5 , 

3 3 15 


(31) 


x = max(n, m) + 1, y = min(n, m) + 1. The derivation of (30) (See Appendix) took into account 
that bosonic terms with non-overlaping sets of indices can be processed simultaneously. 

In the simplified case where = n for each k > K{ er m N s (n, n) = 0(n 5 ) and N g (l,n) = 0(ji 2 ) 
and the single Trotter step overall computational cost is 


0( Y, N 2 N 2 log(N k N l )) + 


, 0<k<l<I<{, 


+oi Y n b N k N?\ + 

V Kf errn <ik^ : l^ : K / 

+o( Y ^logiNMN 2 ). 

\ 0<k<K felln / 


(32) 


. 0 <k<K u 
K {erm <l<I< 


4.0.3 Distinguishable particles 

Distinguishable particles can be distributed to separate classes each with n k = 1 and occupation 
numbers f(p) E {0; 1} with p E S k can be stored using one qubit per spin orbital for each of the 
particles. The creation and annihilation operators can then be represented as 


a t = 

tip ^ _ , 


a_ 


p 

Qjp — G _|_ 


(33) 

(34) 


and we can use either (26) or (32), which are in fact equivalent in this case as n k = 1. 


The equations (26) and (32) demonstrate that the qFCI algorithm can be efficient also for 


systems with bosonic and/or distinguishable particles. In case of the compact boson mapping, the 
smaller number of qubits needed to represent a state of a molecule (as compared with the direct 
boson mapping) is paid by worse computational cost in terms of the number of two-qubit gates 
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needed for a computation. For the first, few-qubit quantum computers, the compact boson mapping 
should be important, however for larger scale quantum computers, the direct boson mapping is of 
much greater importance. 

We would like to note that presented mappings are not restricted to (I)PEA algorithm only. 
They, in fact, can be also used in connection with other methods that were developed to reduce 
qubit and coherence time requirements 21 l 23 l 26 l 55 l, i.e. to adapt the procedure for a present-day or 
near-future quantum technology. 

5 Application to H 2 and HT molecules 

5.1 Computational details 

For the proof-of-principle classical simulations, we have chosen the simplest molecular examples, 
namely the two isotopomers of the hydrogen molecule (Ho, HT). For MOs, we employed the cc- 
pVTZ basis set, while NOs were expanded in a basis of one s, one p, and one d non-contracted 
Gaussians centred on each hydrogen atom with five even-tempered exponents from 9.081045 to 
908.104502 for H and from 27.18608 to 2718.608 for T isotope (total 50 nuclear basis functions for 
each hydrogen atom). The internuclear distance was fixed to R = 0.750746 A. 

The nuclei were for simplicity treated as distinguishable which is the usual procedure that can 
be justified by the fact that exchange interaction between nuclei is negligibly small. 

We worked solely with a compact mapping from subspace of wave functions with zero z- 
component of total electron and nuclear spins and the exponential of a Hamiltonian was simulated 
as an n-qubit gate (similarly as in Refs. 1112 14 15 ). We simulated m = 17 iterations of both IPEA 
A and B with input parameters (see Ref. 14 ) E min = —1.20 a.u. and E max = —1.00 a.u. 

The ground state of both H 2 and HT is dominated by la^la n ilcr n 2 configuration, while the 
excited state u = 1, J = 0 was identified as a state dominated by la^2a n \2a n 2 configuration. 
Between all states in ±50% interval around the experimental transition energy, only this state 
is non-degenerate and symmetric with respect to the exchange of the nuclear coordinates (and 
therefore has even rotation number J and is a nuclear singlet). 

The initial guesses of both states were single determinants (ground state: la^la n \ lrx n 2 , excited 
state: lcr 2 2 cr n i 2 cj ri 2 ). 

However, we must note that for larger polyatomic molecules the identification of several rotation¬ 
less vibrational excited states and the choice of sufficiently accurate initial guesses of eigenvectors 
for IPEA might be much more complicated than it was for the hydrogen molecule. 

5.2 Results 

In Table [lj we show energies and IPEA (A) success probabilities for the ground and excited states 
of H 2 , while Table [2] presents minimal number of repetitions of IPEA A and IPEA B needed to 
achieve a given success probability (0.99, 0.999 999). Tables [3] and [4] give corresponding information 
for the HT molecule. 

The exponential increase of success probabilities p with the number of repetitions r (as supposed 
from the Chernoff bounds 66 ) is demonstrated in Fig. [ 2 ] (for the case of basis consisting of N = 6 
molecular and M = 10 nuclear orbitals) and Fig. [3] ( N = 9, M = 23) by roughly asymptotically 
linear curves of quantity f(r) = — log(l — p) as a function of r. 

For almost all data points, the success probabilities for IPEA A are higher than corresponding 
success probabilities for IPEA B and the slope of f(r) curve in asymptotic region is subsequently 
also higher for IPEA A than for IPEA B. In most cases, we can also see higher values and slopes 
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M 

MOs 

NOs 

E gs 

a.u. 

^es 

a.u. 

A E 

cm -1 

duo 

hgs 

s 

^es 

PA, gs 

PA,es 

1 

4 

(7g 

2(7 7T 

-I.104049 

-1.0809/2 

5064.8 

21.7 

0.9985 

0.9164 

0.8217 

0.8743 

2 

6 

UgU u 

2 otx5 

-1.105483 

-1.082857 

4965.7 

19.3 

0.9973 

0.6737 

0.9645 

0.5760 

3 

8 

2(7 g(7 u 

2(72-1x5 

-1.108131 

-1.085068 

5061.7 

21.6 
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0.7784 

0.8073 

0.7780 

4 

10 

2(7 g2(7 u 

4a2ix 5 

-1.119317 

-1.096120 

5091.1 

22.3 

0.9759 

0.8344 

0.8769 

0.6779 

6 

10 

2 (7g 2 (7 11 7T1 1 


-1.127425 

-1.107187 

4441.8 

6.7 

0.9583 

0.9385 

0.9419 

0.8454 

6 

15 


5a3ix25 

-1.127443 

-1.107373 

4404.8 

5.9 

0.9575 

0.9282 

0.9329 

0.9102 

6 

23 


7a5ir35 

-1.127612 

-1.108051 

4293.3 

3.2 

0.9571 

0.9135 

0.9565 

0.9095 

7 

18 

3(7g2(7 11 7T 11 

6a3-rx3 5 

-1.129852 

-1.109358 

4497.9 

8.1 

0.9404 

0.9041 

0.8113 

0.8056 

7 

23 


7 ct 57 t 3<5 

-1.129917 

-1.110242 

4318.0 

3.8 

0.9464 

0.8900 

0.9059 

0.7301 

9 

10 

3(7g2(7n7rn7Tg 

4(72-7x5 

-1.130100 

-1.109888 

4436.0 

6.6 

0.9419 

0.9053 

0.7718 

0.7387 

9 

18 


6cr3ix35 

-1.130285 

-1.110317 

4382.5 

5.3 

0.9334 

0.8779 

0.7971 

0.8735 

9 

20 


6(74-7x3 5 

-1.130293 

-1.110419 

4361.9 

4.8 

0.9342 

0.8880 

0.9316 

0.8455 

9 

23 


7(75-7x3 5 

-1.130346 

-1.111238 

4193.7 

0.8 

0.9398 

0.8600 

0.8376 

0.7394 
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-1.164025 

-1.145065 

4161.2 
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-1.164033 

-1.145073 

4161.2 








Experiment 


-1.164033 

-1.145073 

4161.2 







Table 1: Ground (E gs ) and excited state ( E es , v = 1, J = 0) NOMO-TRF/FCI energies and IPEA 
(A) success probabilities (PA,gs , PA,es) for the H 2 molecule in the basis consisting of N molecular 
orbitals and M nuclear orbitals. A E denotes transitional energy, 5uo = 100(AE — AE exp ) / AE exp , 
S g s and S es denote square of absolute value of the overlaps between NOMO-TRF/FCI eigenvectors 
and their initial guesses. 


of curves corresponding to ground state when compared to the excited state. This correlates with 
the fact that excited state NOMO-TRF/FCI eigenvector has smaller overlap with its inital guess 
than the ground state. In other words, the excited state has stronger multireference character. 



Figure 2: NOMO-TRF/qFCI success probabilities p (IPEA A - solid line, IPEA B - dashed line) 
for both H 2 and HT molecules in the basis consisting of N = 6 molecular and M = 10 nuclear 
orbitals as a function of the number of repetitions (r). 


5.3 Discussion 

Results presented in the previous section indicate that for all choices of the FCI active space [defined 
by the number of molecular (A) and nuclear (M) orbitals], single determinant initial guesses give 
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Table 2: Minimal number of repetions for IPEA B to achieve success probability at least p = 0.99 
(As,y, 2 ) or at least p = 0.999 999 (Ab^) and the same quantities for IPEA A (AA, y , 2 i NA,yfi), 
where y E {gs(ground state), es(excited state)}. All data corresponds to H 2 molecule. 
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Table 3: Ground (E gs ) and excited state (E es , v = 1, J = 0) NOMO-TRF/FCI energies and IPEA 
(A) success probabilities (pA.gsj PA.es ) for the HT molecule in the basis consisting of N molecular 
orbitals and M nuclear orbitals. A E denotes transitional energy, 5uj = 100(AE — AE exp ) / AE exp , 
S gs and S es denote square of absolute value of the overlaps between NOMO-TRF/FCI eigenvectors 
and their initial guesses. 
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Table 4: Minimal number of repetions for IPEA B to achieve success probability at least p = 0.99 
(-ZVB,y, 2 ) or at least p = 0.999 999 (Ab^g) and the same quantities for IPEA A (IVA,y, 2 , 
where y G {gs(ground state), es(excited state)}. All data corresponds to HT molecule. 



Number of repetitions (r) 


Figure 3: NOMO-TRF/qFCI success probabilities p (IPEA A - solid line, IPEA B - dashed line) 
for both H 2 and HT molecules in the basis consisting of A = 9 molecular and M = 23 nuclear 
orbitals as a function of the number of repetitions (r). 


sufficiently high success probabilities to be amplified by repetitions (p > 0.5). In most cases, success 
probabilities are higher than 0.75. The only exception when the success probability is lower than 
0.6 is the case of the H 2 excited state and N = 2, M = 6. This choice of the FCI active space is 
apparently too small to properly describe the transitional energy anyway. 

Apart from the N = 2, M = 6 active space, results in Tables [2] and [4] indicate that at most 
around 10 repetitions are sufficient to amplify the success probabilities to 0.99, and at most 55 
repetitions for amplification to 0.999999. 

For almost all data points, the success probabilities p for IPEA A are higher than corresponding 
success probabilities for IPEA B and subsequently slopes of f(r) = — log(l — p)) curves in asymp¬ 
totic region are also higher. This could be easily expected as in case of IPEA B, no collapsing of 
the system and improving the overlap between the actual state of the quantum register and the 
exact wave function occurs during iterations which is in contrast to IPEA A. 

In most cases, we can also note higher values and slopes of curves corresponding to ground 
state when compared to the excited state. The only exception is the case of the largest active space 
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used (N = 9 , M = 23) for the HT molecule (Fig. [3]). Lower success probabilities for the excited 
state correlates with the excited state NOMO-TRF/FCI eigenvector having smaller overlap with 
its inital guess than the ground state (stronger multireference character of the excited state) 

In the exceptional case (N = 9, M = 23 for HT), the overlap for the ground state S gs is higher 
than the overlap for the excited state S es , but due to the higher phase reminder 6 (see eqs. 8 and 
9 in 11 ) for the ground state, the success probability of non-repeated (r = 1) IPEA A is slightly 
higher for the excited state than for the ground state. The phase reminder has the same effect for 
repeated IPEA success probabilities in both IPEA A and IPEA B cases. 

6 Conclusions 

In this paper we presented an efficient quantum algorithm for molecular energy computation beyond 
the Born-Oppenheimer approximation. Our approach is based on the quantum full configuration 
interaction method and treats electrons and nuclei on an equal footing, using the nuclear orbital plus 
molecular orbital (NOMO) method. We have presented details of the algorithm and demonstrated 
its performance by simulations on a classical computer. For these simulations we have employed 
relatively small one-particle basis sets and used the compact mapping to keep the number of required 
qubits manageable. 

Two isotopomers of the hydrogen molecule (Ho, HT) were chosen as representative examples 
and calculations of the lowest rotationless vibrational transition energies were simulated. For both 
isotopomers in their ground as well as excited state we have verified that the single-determinant 
initial guess yields high enough success probability to be amplified by repetitions, for both A and 
B version of IPEA. At most 10 repetitions were sufficient to amplify the success probability to 0.99 
and at most 55 repetitions were necessary to achieve 0.999999 success probability. As expected, 
for most data points the success probability of IPEA A was higher than corresponding IPEA B, 
due to the improvement of system wave function overlap due to the measurement in the IPEA A 
procedure. In most cases, the excited state required more repetitions than the ground state, which 
is in agreement with our previous experience on electronic-only calculations^ 14 -! To conclude, the 
qFCI approach has been shown to be viable also for simultaneous treatment of electrons and nuclei 
beyond the Born-Oppenheimer approximation. 


Appendix 

For derivation of the scaling of the compact boson mapping introduced in Section [4j let us consider 
the Hamiltonian parameterised by real-valued integrals V pqrs = V rspq and the most demanding 
4-index term of a single Trotter step 



with p, q, r, and s mutually different. We use the lemma 


exp(ir<I>P) = exp(iT^U DU^) = U exp(iT&D)U\ (36) 

where r and = V pqrs are real numbers, V is a hermitian operator and D = WVU is its diagonal 
form. It is thus sufficient to find eigenvectors and eigenvalues of the operator V = d p d\d s d r + 
aiaid q a p on a Hilbert space corresponding to the quantum register storing occupation numbers 
for spinorbitals p, q, r and s. This space is a direct product of two ./max (p) + 1 = /ma. v(r) + 1 


14 


dimensional and two f max (,Q ) + 1 = fma (s) + 1 dimensional spaces. In the basis characterized by 
boson occupation numbers the operator V has a block-diagonal structure, since 


V\f p , f q J r , fs) = 

= \JUp + 1 ){fq + l)/r/s|/p + 1) fq + 1, fr ~ 1, fs ~ 1) + 

\Jfpfqifr + 1 ){fs + 1)1 fp ~ 1, fq ~ lj fr + 1 5 fs + !)■ 

(37) 

In the simplest case of /max(p) = /max(<z) = ^ there arc 12(n+l—d)+2—d-dimensional blocks 
(d = 1, 2,n+1). In general, the number of diagonal d-dimensional blocks (d = 1, 2,..., min(ni, n 2 )+ 
1) denoted here pa, equals (see the Supplementary Information, Chapter 2.1.) 


Pd = -4,o(N - n 2 \ - l) 2 + 
+2(|ni - n 2 | 2 + 1) + 6z(\ni - n 2 | + z)) 


(38) 


where 2 : = min(m,n 2 ) + 1 - d and m = /max(p) = /max(r) and n 2 = f ma , x (q) = / max (s). De¬ 
composition of the matrix representation of V (V) and then subsequently U into blocks leads to a 
decomposition of U from lemma (36) to a direct sum of unitary operators U l acting on each block 


K 


t> = ©q 


(38) 


i =0 


Let us define 


min(ni,7i2)+l 

M S ,s = dSpd > 


(40) 


d=l 


where s is an auxiliary non-negative integer. Ms iS has a different meaning depending on s value. 
For s = 0, Ms,s equals the total number of subspaces K = Ms ft = nin 2 (rii + n 2 + 1) in the 


decomposition (39). 


The dimension of the quantum register space where operator (35) acts is Mgp = (m + l) 2 (n 2 + 
l) 2 . Ms 1 is in fact the minimal possible size of the quantum register for representing operator 


(35). In the case of qubits, the quantum register dimension will be 2 2 ^ og 2 ( ni+1 )l+ 2 r io S 2 ( n 2 +i)l > 

dimen- 


excess 


(ni + l) 2 (n 2 + l) 2 = Msj ■ Usage of qu-d-its for well chosen d’s may decrease the 
sions usually padded by unit operator blocks and the number C(ni,n 2 ) of classical precomputing 
operations needed for the diagonalization of matrix representation ( |35[ ). 

For s = 2, Mg ]2 = Ng(ni, nf) describes the computational complexity as will be shown in the 
end of this section. 

The computational cost of classical precomputing operations scales as C(ni,n 2 ) = Ms;s- 

Based on the above approach, the exponential (35) can be decomposed into blocks, 


K 


exp(ir4»U) = 


i=0 


Ai = Uiexp(iT$Di)Ui , 


(41) 

(42) 
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where D L is a diagonal gate with only di non-zero elements. Gates A % in the circuit correspond to 
the product (42) and are further described in Fig. [ 5 ] Before and after the sequence of Ai gates is 
applied, the transformation 


I/p; fq, fr, fs) ^ |Si, Ai, A 2 , S 2 ) e->- 

^ |A, Ai, A 2 ,S), 


(43) 


and its inverse have to be applied as shown in Fig. [4] (The transformation (43) is realized by the 
subcircuit in the dashed box (Fig. [ 4 J the gate W). 

The non-negative integers f p , f q , f r and f s are occupation numbers, their upper-bounds are 
f p , f r < n\ and f q , f s < n 2 (and corresponding register sizes in qubits Q 1 = (log 2 (ni + 1)] and 


Qi = flog 2 (n *2 + 1)])- In formula (43) the first transformation produces the sum and difference of 
occupation number pairs (/ p , f q ) and (/ r , f s ), 


Si = fq + fp, 

(44) 

> 

II 

1 

"P 

(45) 

S 2 = fs + fr, 

(46) 

< 

1 

II 

< 

(47) 


and the second transformation in (43) produces sum and difference of the first and third registers 
denoted £1 and £2 respectively, 


£ — £2 + £1 


(48) 


A = £ 2 —£ 1 . (49) 

The whole Trotter step (35) is represented by a quantum circuit model in Fig. [4] and the par¬ 
ticular implementation of the ASG gate is discussed in detail in Chapter 2.3 of the Supplementary 
Information. Note that the ASG gate can be realized by O(Q) elementary gates and (depending 
on a particular realization) with 0 to Q working qubits. The existence of the ASG gate is obvious 
from the fact that from the combination of the sum and difference of the two input integers, the 
input integers can be unambiqously deduced. The ASG gate outputs for the ancilla registers | a*) 
different from zeros are irrelevant. 
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I°l) 
I fp) 

\fg) 

l a s) 

I fr) 
\fs) 
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Figure 4: The quan tum circuit implementing a single Trotter step (boson-boson interaction, see Eq. 
(35), (36) and (39)). The Ai Q + 2 qubit gates are controlled by ^-specific Q + 1, Q + 2 and Q + 1 
bit sequences, quantum circuit corresponding to them is presented in Fig. [5] with one additional 
working qubit increasing the number of qubits Ai operates on to Q + 3. The gate labeled ASG 
is the Adder-Subtractor Gate realizing the simultaneous addition (sum) and subtraction of two 
binary represented integers. The different indices for ASG gates corresponds to different ordering 
of outputs and inputs. The \Q\ — Q2I + 2-qubit kets |ai) and I02) are initialy set to zero (|0)) and 
correspond to the neccessity to pad the input integer with less bits with |Qi — Q 2 I zeros to match 
the number of bits of the greater input integer prior the particular logical operation, one carry bit 
is needed for the output integer corresponding to the sum and one extra bit (sign bit) is needed for 
the output integer corresponding to the difference. Similarly, the ket ( 03 ) is a 2-qubit initialy set 
to zero (| 00 )) to provide the most significant bit for the sum £ and the sign bit for the difference 
A. The ket ( 04 ) is a R = 2|Qi — Q 2 I + 6 register for the recovery of all working qubits. 
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l(A)i) 

l(A) 2 > 

l(A)a) 


(A)[log 2 rfil + 2 ) 


(A)[io g2 41+3) 


|(A)q +2 ) 

|(A)q +3 ) 



Figure 5: The quantum circuit representing the action of the gate (which occures in Fig. [ 4 ] and 
equations (41) and (42)) on its target quantum register A padded by one ancilla qubit (A)q + 3 
needed for ADD gate operation. The Q+3 qubit gate ADD (a) adds a constant integer a = ±Ao,i 
to the input quantum register producing the output quantum register. The ADD gate is based 
on the generalized <f> ADD gate (operating in 0(1) time, see the Fig. 4 on the page 654 of 67 '' 


inserted between the forward and backward Quantum Fourier Transforms (operating in O(QlogQ) 
time each). The value Ao,* corresponding to the z-th subspace (in the decomposition (p9|)) is the 
smallest value of A corresponding to any vector in that subspace and is to be precomputed classicaly 
for each i in {1, 2,..., K}. The primes correspond to the fact that the operators now act on different 


subspaces. In the special case of one dimensional subspaces (d, = 1), U-, exp(zr<f>l3') and gates 
have 0 qubit target register and are therefore equal to multiply controlled phase-shifts. Since the 
corresponding phase is 1 = exp(O) for all of them, they can be omitted from the quantum circuit. 
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The transformation (43) needs 2\Q± — Q2I + 4 ancilliary qubits for the first part and another 2 
ancilliary qubits for the second part. Therefore the first and the last quantum registers (storing 
the values of A and £ respectively) have a size max(Qi, Q 2 ) + 2 qubits while the middle quantum 
registers (storing the values of Ai and A 2 ) of max(Q\, Q 2 ) + 1 qubits. In the following text, 
Q = max((5i, Qi)- The number of single qubit gates and CNOTs for the transformation ( |43| ) 
scales as 0(Q ) if the algorithm presented in 68 is used (with another Q + 1 working qubits) or 
as 0(QlogQ ) (but with no need for further working qubits) if the algorithm exploiting QFT®,®1 
is employed. Each Ui quantum gate acts on the target register which is subregister of the first 
register in the last part of (43), storing the value of A and is multiply controlled by other qubits 
representing the ket |A,Ai,A 2 ,£). Multiply controlled quantum gates can be decomposed into 
the bare quantum gate acting on the target register and either O(Q) (1-qubit and CNOT) gate 
cost with using 3Q + 2 working ancilliary qubits or (D(Q 2 ) gate cost without any working ancilliary 
qubits (for the algorithm see- 1 on pages 183 and 193 respectively, the variant with ancilliary qubits 
is also mentioned in' 1 ). The action of Ai gate from Fig. [ 4 ] and equations ( |4l| ) and (42) is described 


in the quantum circuit in Fig. |5j The bare quantum gate from Fig. |5| U-, acting on the target 
register (2 r io sU)l -dimensional subspace) can be decomposed (via Quantum Shanon Decomposition 
(QSD),^ 2 ) into 0(d 2 ) elementary quantum gates. Neglecting the contribution of gates acting on 
the controlling register which scales as 0(MsflQ 2 ) = 0(x 2 ylog 2 (x + 1)) in the worst case, the 
formula for A gates (n, m) (31) is derived as Mg t 2 - The diagonal bare quantum gate exp(zr<I>.D') can 
be decomposed into the 1-qubit gates and CNOTs at most with the same effort as U[ as further 
discussed in the Supplementary Information (Chapter 2.4.). It is important to note, that the 
classical pre-processing - diagonalization of V and QSD of the corresponding Ui operators needs 
to be done just once (for e.g. p = 1, q = 2,r = 3, s = 4) before the quantum algorithm is started, 
then for each elementary Trotter term (35) the seqence of quantum gates differs just by addressing 
different quantum registers (no longer p = l,g = 2,r = 3,s = 4) and by different value of 4> = V pqrs 
as term in phase-parameter in exp(ir&D) diagonal operator from lemma (36). 
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